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Abstract 



A direct method for the computation of polynomial conservation laws of polynomial systems 
of nonlinear partial differential equations (PDEs) in multi-dimensions is presented. The 
method avoids advanced differential-geometric tools. Instead, it is solely based on calculus, 
variational calculus, and linear algebra. 

Densities are constructed as linear combinations of scaling homogeneous terms with un- 
determined coefficients. The variational derivative (Euler operator) is used to compute the 
undetermined coefficients. The homotopy operator is used to compute the fluxes. 

The method is illustrated with nonlinear PDEs describing wave phenomena in fluid dy- 
namics, plasma physics, and quantum physics. For PDEs with parameters, the method 
determines the conditions on the parameters so that a sequence of conserved densities might 
exist. The existence of a large number of conservation laws is a predictor for complete integra- 
bility. The method is algorithmic, applicable to a variety of PDEs, and can be implemented 
in computer algebra systems such as Mathematica, Maple, and REDUCE. 

1 Introduction 

Nonlinear partial differential equations (PDEs) that admit conservation laws arise in many 
disciplines of the applied sciences including physical chemistry, fluid mechanics, particle 
and quantum physics, plasma physics, elasticity, gas dynamics, electromagnetism, magneto- 
hydro-dynamics, nonlinear optics, and the bio-sciences. Conservation laws are fundamental 
laws of physics. They maintain that a certain quantity, e.g. momentum, mass (matter), 
electric charge, or energy, will not change with time during physical processes. Often the 
PDE itself is a conservation law, e.g. the continuity equation relating charge to current. 

As shown in [8] and the articles in this issue, computer algebra systems (CAS) like Math- 
ematical Maple, and REDUCE, are useful to tackle computational problems in chemistry. 
Finding closed-form conservation laws of nonlinear PDEs is a nice example. Using CAS 
interactively, we could make a judicious guess (ansatz) and find a few simple densities and 
fluxes. Yet, that approach is fruitless for complicated systems with nontrivial conservation 
laws. Furthermore, completely integrable PDEs [1, 2] admit infinitely many independent 
conservation laws. Computing them is a challenging task. It involves tedious computations 
which are prone to error if done with pen and paper. The most famous example is the 
Korteweg-de Vries (KdV) equation from soliton theory [1, 2] which describes water waves 
in shallow water, ion-acoustic waves in plasmas, etc. Our earlier work [3, 16] dealt primarily 
with the symbolic computation of conservation laws of completely integrable PDEs in (1 + 1) 
dimensions (with independent variables x and t). In this paper we present a symbolic method 
that covers PDEs in mult i- dimensions (e.g. x, y, z, and t), irrespective of their complete inte- 
grability. As before, our approach relies on the concept of dilation (scaling) invariance which 
limits it to polynomial conserved densities of polynomial PDEs in evolution form. 

There are many reasons to compute conserved densities and fluxes of PDEs explicitly. 
Invariants often lead to new discoveries as was the case in soliton theory. We may want to 
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verify if conserved quantities of physical importance (e.g. momentum, energy, Hamiltonians, 
entropy, density, charge) are intact after constitutive relations have been added to close a 
system. For PDEs with arbitrary parameters we may wish to compute conditions on the 
parameters so that the model admits conserved quantities. Conserved densities also facilitate 
the study of qualitative properties of PDEs, such as bi- or tri-Hamiltonian structures. They 
often guide the choice of solution methods or reveal the nature of special solutions. For 
example, an infinite sequence of conserved densities assures complete integrability [1] of the 
PDE, i.e. solvability by the Inverse Scattering Transform [1] and the existence of solitons [2]. 

Conserved densities aid in the design of numerical solvers for PDEs [27]. Indeed, semi- 
discretizations that conserve discrete conserved quantities lead to stable numerical schemes 
(i.e. free of nonlinear instabilities and blowup). While solving differential-difference equations 
(DDEs), which arise in nonlinear lattices and as semi-discretizations of PDEs, one should 
check that their conserved quantities indeed remain unchanged. Capitalizing on the analogy 
between PDEs and DDEs, the techniques presented in this paper have been adapted to 
DDEs [12, 17, 19, 20, 21] and fully discretized lattices [14, 15]. 

There are various methods (see [19]) to compute conservation laws of nonlinear PDEs. A 
common approach relies on the link between conservation laws and symmetries as stated in 
Noether's theorem [6, 22, 26]. However, the computation of generalized (variational) sym- 
metries, though algorithmic, is as daunting a task as the direct computation of conservation 
laws. Nonetheless, we draw the reader's attention to DE_APPLS, a package for constructing 
conservation laws from symmetries available within Vessiot [7], a general purpose suite of 
Maple packages for computations on jet spaces. Other methods circumvent the existence of 
a variational principle [4, 5, 29, 30] but still rely on the symbolic solution of a determining 
system of PDEs. Despite their power, only a few of the above methods have been fully 
implemented in CAS (see [16, 19, 29]). 

We purposely avoid Noether's theorem, pre-knowledge of symmetries, and a Lagrangian 
formulation. Neither do we use differential forms or advanced differential- geometric tools. 
Instead, we present and implement our tools in the language of calculus, linear algebra, 
and variational calculus. Our down-to-earth calculus formulas are transparent, easy to use 
by scientists and engineers, and are readily adaptable to nonlinear DDEs (not covered in 
Vessiot) . 

To design a reliable algorithm, we must compute both densities and fluxes. For the latter, 
we need to invert the divergence operator which requires the integration (by parts) of an 
expression involving arbitrary functions. That is where the homotopy operator comes into 
play. Indeed, the homotopy operator reduces that problem to a standard one-dimensional 
integration with respect to a single auxiliary parameter. One of the first 1 uses of the ho- 
motopy operator in the context of conservation laws can be found in [5]. The homotopy 
operator is a universal, yet little known, tool that can be applied to many problems in which 
integration by parts in multi- variables is needed. A literature search revealed that homotopy 
operators are used in integrability testing and inversion problems involving PDEs, DDEs, 
lattices, and beyond [19]. Assuming the reader is unfamiliar with homotopy operators, we 

1 We refer the reader to [26, p. 374] for a brief history of the homotopy operator. 
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"demystify" the homotopy formulas and make them ready for use in nonlinear sciences. 

The particular application in this paper is computation of conservation laws for which 
our algorithm proceeds as follows: build a candidate density as a linear combination (with 
undetermined coefficients) of "building blocks" that are homogeneous under the scaling sym- 
metry of the PDE. If no such symmetry exists, construct one by introducing parameters with 
scaling. Subsequently, use the conservation equation and the variational derivative to derive 
a linear algebraic system for the undetermined coefficients. After the system is solved, use 
the homotopy operator to compute the flux. Implementations for (1 + l)-dimensional PDEs 
[16, 17] in Mathematica and Maple can be downloaded from [10, 18]. Mathematica code that 
automates the computations for PDEs in multi-dimensions is being designed. 

This paper is organized as follows. We establish some connections with vector calcu- 
lus in Section 2. In Section 3 we list the nonlinear PDEs that will be used throughout 
the paper: the Korteweg-de Vries and Boussinesq equations from soliton theory [1, 2], the 
Landau-Lifshitz equation for Heisenberg's ferromagnet [13] and a system of shallow water 
wave equations for ocean waves [1 1] . Sections 4 and 5 cover the dilation invariance and con- 
servation laws of those four examples. In Section 6, we introduce and apply the tools from 
the calculus of variations. Formulas for the variational derivative in multi-dimensions are 
given in Section 6.1, where we apply the Euler operator for testing the exactness of various 
expressions. Removal of divergence-equivalent terms with the Euler operator is discussed 
in Section 6.2. The higher Euler operators and homotopy operators in ID and 2D are in 
Sections 6.3 and 6.4. In the latter section, we apply the homotopy operator to integrate 
by parts and to invert of the divergence operator. In Section 7 we present our three-step 
algorithm to compute conservation laws and apply it to the KdV equation (Section 7.1), 
the Boussinesq equation (Section 7.2), the Landau-Lifshitz equation (Section 7.3), and the 
shallow water wave equations (Section 7.4). Conclusions are drawn in Section 8. 

2 Connections with Vector Calculus 

We address a few issues in multivariate calculus which complement the material in [9]: (i) 
To determine whether or not a vector field F is conservative, i.e. F = V/ for some scalar 
field /, one must verify that F is irrotational or curl free, that is V x F = 0. The cross (x) 
denotes the Euclidean cross product. The field / can be computed via standard integrations 
[24, p. 518, 522]. 

(ii) To test if F is the curl of some vector field G, one must check that F is incompressible 
or divergence free, i.e. V • F = 0. The dot (■) denotes the Euclidean inner product. The 
components of G result from solving a coupled system of first-order PDEs [24, p. 526]. 

(iii) To verify whether or not a scalar field / is the divergence of some vector function F, no 
theorem from vector calculus comes to the rescue. Furthermore, the computation of F such 
that / = V ■ F is a nontrivial matter which requires the use of the homotopy operator or 
Hodge decomposition [9]. In single variable calculus, it amounts to computing the primitive 
F = ffdx. 
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In multivariate calculus, all scalar fields /, including the components Fj of vector fields 
F = (F ± , F 2 , F 3 ), are functions of the independent variables (x, y, z). In differential geometry 
the above issues are addressed in much greater generality. The functions / and Fi can now 
depend on arbitrary functions u(x,y, z),v(x,y, z), etc. and their mixed derivatives (up to 
a fixed order) with respect to (x,y,z). Such functions are called differential functions [26]. 
As one might expect, carrying out the gradient-, curl-, or divergence-test requires advanced 
algebraic machinery. For example, to test whether or not / = V • F requires the use of 
the variational derivative (Euler operator) in 2D or 3D as we will show in Section 6.1. The 
actual computation of F requires integration by parts. That is where the higher Euler and 
homotopy operators enter the picture (see Sections 6.3 and 6.4). 

At the moment, no major CAS have reliable routines for integrating expressions involving 
arbitrary functions and their derivatives. As far as we know, CAS offer no functions to test 
if a differential function is a divergence. Routines to symbolically invert the total divergence 
are certainly lacking. In Section 6 we will present these tools and apply them in Section 7. 



3 Examples of Nonlinear PDEs 

Definition: We consider a nonlinear system of evolution equations in (3 + 1) dimensions, 
u t = G(u, 

Uxj Uy, U 2 , U 2 xj U-2y, U.22, U. X y, W X zi Hyzi ■ ■ ■ ), (1) 

where x= (x, y, z) and t are space and time variables. The vector u(x, y, z, f) has N compo- 
nents U{. In the examples we denote the components of u by u,v,w, etc.. Throughout the 
paper we use the subscript notation for partial derivatives, 

du d 2 u d 5 u d 2 u 

Ut= dt' U2x = Uxx= dx 1, u ^y = u ™yyy = d^dyl> u *y = etc - (2) 

We assume that G is smooth [9] and does not explicitly depend on x and t. There are no 
restrictions on the number of components, order, and degree of nonlinearity of G. 

We will predominantly work with polynomial systems, although systems involving one 
transcendental nonlinearity can also be handled [3, 19]. If parameters are present in (1), 
they will be denoted by lower-case Greek letters. Throughout this paper we work with the 
following four prototypical PDEs: 

Example 1: The ubiquitous Korteweg-de Vries (KdV) equation [1, 2] 

u t + uu x + u 3x = 0, (3) 

for u(x,t) describes unidirectional shallow water waves and ion-acoustic waves in plasmas. 
Example 2: The wave equation, 

u a ~ u 2x + 3uu 2x + 3u 2 x + au Ax = 0, (4) 
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for u(x,t) with real parameter a, was proposed by Boussinesq to describe surface waves in 
shallow water [1]. For what follows, we rewrite (4) as a system of evolution equations, 

u t + v x = 0, v t + u x - 3uu x - au 3x = 0, (5) 

where v(x,t) is an auxiliary dependent variable. 
Example 3: The classical Landau-Lifshitz (LL) equation, 

S t = S x AS + S x DS, (6) 

without external magnetic field, models nonlinear spin waves in a continuous Heisenberg 
ferromagnet [13, 28]. The spin vector S(x, t) has real components (u(x,t),v(x,t),w(x,t)); 
A = V 2 is the Laplacian, and D = diag(a, /3, 7) is a diagonal matrix with real coupling 
constants a, (3, and 7 along the downward diagonal. Splitting (6) into components we get 

Ut = vw 2x - wv 2x + (7 - P)vw, 

v t = wu 2x - uw 2x + {a - ~f)uw, (7) 
w t = uv 2x - vu 2x + ((3 - a)uv. 

The second and third equations can be obtained from the first by cyclic permutations 

(u — > v, v — > w, w — > u) and (a — > f3, (5 — > 7, 7 — > a). (8) 

We take advantage of the cyclic nature of (7) in the computations in Section 7.3. 
Example 4: The (2+l)-dimensional shallow- water wave (SWW) equations, 

u t + (u-V)u + 2Q x u = - V(/i0) + \hV6, 

ht + V-(/m) = 0, t + u-(V0)=O, (9) 

describe waves in the ocean using layered models [11]. Vectors u = u(x, y,t)i + v(x, y,t)j 
and Q = f2k are the fluid and angular velocities. V = + is the gradient operator and 
i, j, and k are unit vectors along the x, y, and z-axes. 8(x,y,t) is the horizontally varying 
potential temperature field and h(x, y, t) is the layer depth. System (9) can be written as 

u t + uu x + vuy — 2Vtv + \h6 x + 6h x = 0, v t + uv x + vv y + 2Qu + \h6 y + 6h y = 0, 

h t + hu x + uh x + hv y + vh y = 0, 9 t + u6 x + v6 y — 0. (10) 

4 Key Concept: Dilation or Scaling Invariance 

Definition: System (1) is dilation invariant if it does not change under a scaling symmetry. 
Example: The KdV equation (3) is dilation invariant under the scaling symmetry 

(x, t, u) -> (X^x, \~H, \ 2 u), (11) 



5 



where A is an arbitrary parameter. Indeed, apply the chain rule to verify that A 5 factors out. 

Definition: The weighft, W, of a variable is the exponent in X p which multiplies the variable. 

Example: We will always replace x by \~ 1 x. Thus, W(x) = —1 or W(d/dx) = 1. In view 
of (11), W(d/dt) = 3 and W(u) = 2 for the KdV equation (3). 

Definition: The rank of a monomial is defined as the total weight of the monomial. An 
expression is uniform in rank if its monomial terms have equal rank. 

Example: All monomials in (3) have rank 5. Thus, (3) is uniform in rank with respect to 
(11)- 

Weights of dependent variables and weights of d/dx,d/dy, etc. are assumed to be non- 
negative and rational. Ranks must be positive integer or rational numbers. The ranks of the 
equations in (1) may differ from each other. 

Conversely, requiring uniformity in rank for each equation in (1) we can compute the 
weights, and thus the scaling symmetry, with linear algebra. 

Example: Indeed, for the KdV equation (3) we have 

W{u) + W(d/dt) = 2W(u) + 1 = W(u) + 3, (12) 
where we used W(d/dx) = 1. Solving (12) yields 

W{u) = 2, W(d/dt) = 3. (13) 

This means that x scales with A _1 ,i with A~ 3 , and u with A 2 , as claimed in (11). 

Dilation symmetries, which are special Lie-point symmetries [26], are common to many 
nonlinear PDEs. However, non-uniform PDEs can be made uniform by giving appropriate 
weights to parameters that appear in the PDEs, or, if necessary, by extending the set of 
dependent variables with auxiliary parameters with appropriate weights. Upon completion 
of the computations we can set the auxiliary parameters equal to 1. 

Example: The Boussinesq system (5) is not uniform in rank because the terms u x and au 3x 
lead to a contradiction in the weight equations. To circumvent the problem we introduce an 
auxiliary parameter (3 with (unknown) weight, and replace (5) by 

u t + v x = 0, 

v t + (3u x - 3uu x - au 3x = 0. (14) 
Requiring uniformity in rank, we obtain (after some algebra) 

W(u) = 2, W(v) = 3, W(0)=2, W(^) = 2. (15) 

Therefore, (14) is invariant under the scaling symmetry 

(x, t, u, v, (5) -> (A^rr, \-\ \ 2 u, X s v, \ 2 (3). (16) 
2 The weights are the remnants of physical units after non-dimensionalization of the PDE. 
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Definition: System (1) is called multi-uniform in rank if it admits more than one scaling 
symmetry (none of which result from introducing unnecessary parameters with weights). 

Example: The LL system (7) is not uniform in rank unless we allow the parameters a, (3, 
and 7 to have weights. Doing so, uniformity in rank requires 

d 

W(u)+W(—) = W(v) +W(w) +2 = W(j)+W{v)+W{w) = W{(3) + W(v) + W(w) , 
W(v) + W(—) = W(u) + W(w)+2 = W(a)+W(u) + W(w) = W(-f) + W(u) + W(w), (17) 

C6 

d 

W(w)+W(— ) = W(u)+W(v)+2 = W{/3)+W(u)+W(v)=W(a)+W(u)+W(v). 

Solving these relations gives W(u) = W(v) = W(w), W(a) = W{(3) = W( 7 ) = 2, W(§- t ) = 
W{u) + 2, where W{u) is arbitrary. The LL equation is thus multi-uniform and invariant 
under the following class of scaling symmetries, 

(x, t, u, v, w, a, (3, 7) -> (X^x, A" (a+2) t, X a u, X a v, X a w, X 2 a, X 2 (3, A 2 7 ), (18) 

where W{u) = a is an arbitrary non-negative integer or rational number. For example, if 
we take a = 1 then W{u) = W(v) = W(w) = 1, W(a) = W{0) = W(i) = 2, W(f t ) = 
3 and (7) is invariant under (x, t, u, v, w, a, (5, 7) — > (A _1 a;, X~H, Xu, Xv , Xw, X 2 a, X 2 j3, X j). 
Taking different choices for a < W(a) = 2 will prove advantageous for the computations in 
Section 7.3. 

Example: The SWW equations (10) are not uniform in rank unless we give a weight to Q. 
Indeed, uniformity in rank for (10) requires, after some algebra, that 

W(d/dt) = W{tt), W(d/dy) = W(d/dx) = 1, W(u) = W(v) = W{tt) - 1, 
W(6)=2W(Sl)-W(h)-2, (19) 

where W(h) and are arbitrary. Hence, the SWW system is multi-uniform and invariant 

under the class of scaling symmetries, 

(x, y, t, u, v, h, 6, Q) -> (X~ 1 x, X~ 1 y, X~ b t, X^u, X^v, X a h, X 2b ~ a ~ 2 9, X b Q). (20) 

where W(h) = a and W(Q) = b. Various choices for a and b will aid the computations in 
Section 7.4. 



5 Conservation Laws 

Definition: A conservation law [26] in differential form is a PDE of the form, 

D t p + DivJ = 0, (21) 

which is satisfied on solutions of (1). The total time derivative D t and total divergence 
Div are defined below. The scalar differential function p is called the conserved density; 
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the vector differential function J is the associated flux. In electromagnetism, (21) is the 
continuity equation relating charge density p to current J. 

In general, p and J are functions of x, t, u, and the derivatives of u with respect to the 
components of x. In this paper, we will assume that densities and fluxes (i) are local, which 
means free of integral terms; (ii) polynomial in u and its derivatives; and (iii) do not explicitly 
depend on x and t. So, p(u, u^,, u y , u z , u 2x , . . .) and J(u, u x , u y , u z , u 2x , . . .). 

Definition: There is a close relationship between conservation laws and constants of motion. 
Indeed, integration of (21) over all space yields the integral form of a conservation law, 

/'OO poo roo 
/ / p dxdydz = constant, (22) 
-oo J—oo J — OO 

provided that J vanishes at infinity. As in mechanics, the P's are called constants of motion. 

The total divergence Divis computed as Div J = (D x , D^, D 2 )-(Ji, J 2 , J 3 ) = D x Ji+D y J 2 + 
D Z J 3 . In the ID case, with one spatial variable (x), conservation equation (21) reduces to 

D t p + D x J = 0, (23) 

where both density p and flux J are scalar differential functions. 

The conservation laws (21) and (23) involve total derivatives D t and D^,. In the ID case, 

D ^t+£^ D ' W - (24) 

where n is the order of u in p. Upon replacement of u t , u tx , etc. from u t = G, we get 

B t p=^ + p(uy[G], (25) 
where p(u)'[G\ is the Frechet derivative of p in the direction of G. Similarly, 

dJ m dJ 

D - J= & + £s^"" +i »- (26) 

where m is the order of u of J. 

Example: Assume that the KdV equation (3) has a density of the form 

p = du 3 + c 2 u 2 x , (27) 

where c\ and c 2 are constants. Since p does not explicitly depend on time and is of order 
n — 1 in u, application of (24) gives 

= 3ciu 2 u t + 2c 2 u x u tx = -3cim 2 (mm :e + u 3x ) - 2c 2 u x (uu x + u 3x ) x 
= -?,ciu 2 (uu x + u 3x ) - 2c 2 u x (u 2 x + uu 2x + u 4x ) 

= -(3c 1 u 3 u x + 3ciu 2 u 3x + 2c 2 u x + 2c 2 uu x u 2x + 2c 2 u x u 4:X ), (28) 
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where D° = I is the identity operator and where we used (3) to replace u t and u tx . The 
generalizations of (24) and (26) to multiple dependent variables is straightforward [26]. 

Example: Taking u(x, t) = (u(x,t),v(x,t),w(x,t)), 

dJ mi dJ 1112 dJ m3 dJ 

DxJ = ^ + E oT^Hilx + E + E (30) 

OX k=Q OUkx fc=0 CVkx k=Q OWkx 



where rii,n 2 and n 3 are the highest orders of u, v and w in p, and m 1? m 2 and m 3 are the 
highest orders of u, v , w in J. 

Example: Assume that the Boussinesq system (14) has a density of the form 

p = cx(3 2 u + c 2 /3u 2 + c 3 u 3 + C4V 2 + c 5 u x v + c e ul. (31) 
where C\ through c 6 are constants. Since n\ — \ and n 2 = 0, application of (29) gives 

= (ci/3 2 + 2c 2 (3u + 3c 3 u 2 )u t + (c 5 v + 2c 6 u x )u tx + {2c A v + c 5 u x )v t 
= -((c 1 (3 2 + 2c 2 [3u + 3c 3 u 2 )v x + (c 5 v + 2c 6 u x )v 2x + (2c A v + c 5 u x )(^^ 

where we replaced u t ,u tx , and v t from (14). 

Remark: The flux J in (21) is not uniquely defined. In ID we use (23) and the flux is only 
determined up to an arbitrary constant. In 2D, the flux is only determined up to a divergence- 
free vector K = (Ki, K 2 ) = (D,^, — D x <f>), where <fi is an arbitrary scalar differential function. 
In 3D, the flux can only be determined up to a curl term. Indeed, if (p, J) is a valid density- 
flux pair, so is (p, J + V x K) for any arbitrary vector differential function K = (Ki, K 2 , K 3 ). 
Recall that V x K = (D y K 3 - B Z K 2 , B z Ki - D X K 3 , D X K 2 - DyK^. 

Example: The Korteweg-de Vries equation (3) is known [1, 2] to have infinitely many 
polynomial conservation laws. The first three density-flux pairs are 



(1) = u, JW = \u 2 + u 2x , 

1 , ,™ 1 , 1 



P (2) = ^i 2 , J {2) = V - tu 2 x + uu 2x , (33) 

1 1 

p^ = -u 3 — u 2 x) = -M 4 — 2uu 2 x + u 2 u 2x + u\ x — 2u x u 3x . 

The first two express conservation of momentum and energy, respectively. They are easy 
to compute by hand. The third one, less obvious and requiring more work, corresponds to 
Boussinesq's moment of instability [25]. Observe that the above densities are uniform of 
ranks 2,4, and 6. The fluxes are also uniform of ranks 4, 6, and 8. 
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The computations get quickly out of hand for higher ranks. For example, for rank 12 



P {6) = - 10w 3 w 2 - hu x + 18u 2 u 2x 2 + ™u 2x 3 - ^uu 3x 2 + yM 4 x 2 (34) 

and is a scaling homogeneous polynomial with 20 terms of rank 14 (not shown). 
In general, if in (23) rank(p) = R then rank(J) = R + W(d/dt) - 1. All the pieces in (21) 
are also uniform in rank. This comes as no surprise since the conservation law (21) holds on 
solutions of (1), hence it 'inherits' the dilation symmetry of (1). 

Example: The Boussinesq equation (4) also admits infinitely many conservation laws and 
is completely integrable [1, 2]. The first four density-flux pairs [3] for (14) are 



P (1) = u, 

= v , =p u - -u 2 - au 2x , (35) 

2 

= uv, J*- 3 -* = -f3u 2 — u 3 H — v 2 H — au 2 — auu 2x , 

2 2 2 

p^ = f3u 2 — u 3 + v 2 + au 2 , = 2j3uv — 3u 2 v — 2au 2x v + 2au x v x . 

The densities are of ranks 2, 3, 5 and 6, respectively. The corresponding fluxes are of one 
rank higher. After setting (3 = 1 we obtain the conserved quantities of (5) even though 
initially this system was not uniform in rank. 

Example: We assume that a, (3, and 7 in (6) are nonzero. Cases with vanishing parameters 
must be investigated separately. The first six density-flux pairs [13, 28] are 





= u, 




= v x w — vw x , 


W 


= 7 ^ a), 


p (2) 


= v, 




= uw x — u x w, 


(a 


= 7^), 


p (3) 


= w, 




J (3) = u x v - uv x , 


(a 




p (4) 


= u 2 + v 2 


+ w 2 , 


J< 4 > = 0, 






p (5) 


= (u 2 + v 


2 + w 2 ) 2 , 


= 0, 






p (6) 


2 1 2 

= u x + v x 


+ ™l + (7 - a) 








J(6) 


= l({yw x 


- v x w)u 2x + (u. 


x w - uw x )v 2x + (uv x 


— 11 


xV)w 2x 




+ 09- 


- j)u x vw + (7 - 


a)uv x w + (a — (3)uvw x 


I- 



(36) 



Note that the second and third density-flux pairs follow from the first pair via the cyclic 
permutations in (8). If we select W(u) = W(v) = W(w) — a — |, then the first three 
densities have rank | and p( 4 \p( 5 \ and p^ have ranks |, 1, §, respectively. 

The physics of (6) demand that the magnitude S — ||S|| of the spin field is constant in 
time. Indeed, using (7) we readily verify that uu t + vv t + ww t = 0. Hence, any different iable 
function of S is also conserved in time. This is apparent in p^ = ||S|| 2 and p^ = ||S|| 4 , 
which are conserved (in time) even for the fully anisotropic case where a^/3^7^«. 
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Obviously, linear combinations of conserved densities are also conserved. Hence, p*- 6 -* 
\Sx\\ 2 + (7 — cx)u 2 + (7 — (3)v 2 can be replaced by 



= ^(P (6) " 7P (4) ) = \{u 2 x + vl + wl - an 2 - (3v 2 - W)- (37) 



Now, u\ = D x (uu x ) — uu2x, etc., and densities are divergence-equivalent (see Section 6.2) if 
they differ by a total derivative with respect to x. So, (37) is equivalent with 

p (6) = -^(mm 2x + vv 2x + ww 2x + au 2 + (3v 2 + 7W 2 ), (38) 
which can be compactly written as 

p(6) = _I(S- AS + S-DS), (39) 
2 

where A is the Laplacian and D = diag(a, f3, 7). So, is the Hamiltonian density [28] in 

H = ~ J (S • AS + S • DS) dx 

= I J (\\S X \\ 2 + - a)u 2 + ^ - (3)v 2 ) dx. (40) 

The Hamiltonian H expresses that the total energy is constant in time since dH/dt = 0. 

Example: The SWW equations (10) admit a Hamiltonian formulation [11] and infinitely 
many conservation laws. Yet, that is still insufficient to guarantee complete integrability of 
PDEs in (2 + l)-dimensions. The first few conserved densities and fluxes for (10) are 

P w = h, j« ( "'! ) . P ( 2 ) = he, j( 2 ) 1 " ,lf) ■ 



p(?) = h 2 , j( 3 ) 



vh J ' ' \vh9 

uh6 2 \ 
vh6 2 J ' 



(4) _ (u 2 , v 2 )h , ^ j(4) _ ( « 3 /l + UV 2 h + 2uh 2 9 \ ( 41 ) 



P ( 5 ) = (2fi - Uy + 7^)0, 



(5) _ 1 / \2Vtu9 - AuUyO + 6uu x + 2vv y 9 + m 2 ^ + v 2 y - h99 y + /i,,6> 2 \ 
" 6 ^ I2ttv9 + Avv x 9 - 6vu y 9 - 2uu x 9 - u 2 9 x -v 2 9 x + h99 x - h x 9 2 ) ' 

As shown in [11], system (10) has conserved densities p = hf(9) and p = (20 — u y + v x )g(9), 
where f{9) and g{9) are arbitrary functions. Such densities are the integrands of the Casimirs 
of the Poisson bracket associated with (10). The algorithm presented in Section 7.4 only 
computes densities of the form p = h9 k and p = (2Q — u y + v x )9 l , where k and / are positive 
integers. 
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6 Tools from the Calculus of Variations 



In this section we introduce the variational derivative (Euler operator), the higher Euler 
operators (also called Lie-Euler operators) from the calculus of variations, and the homotopy 
operator from homological algebra and variational bi-complexes [6]. These tools will be 
applied in Section 7. 

6.1 Variational Derivative (Euler Operator) 

Definition: A scalar differential function / is a divergence if and only if there exists a 
vector differential function F such that / = DivF. In ID, we say that a differential function 
/ is exact? if and only if there exists a scalar differential function F such that / = D X F. 
Obviously, F = D~ 1 (/) = J f dx is then the primitive (or integral) of /. 

Example: Taking / = — D t p from (28), that is 

/ = 3ciu 3 u x + 2c 2 u 3 x + 2c 2 uu x u 2x + 3ciU 2 u 3x + 2c 2 u x u 4x , (42) 

where 4 u(x; t), we will compute c\ and c 2 so that / is exact. 
Upon repeated integration by parts (by hand), we get 

F — J f dx = ^ci-u 4 + (c 2 - 3ci)uu 2 x + 3c 1 u 2 u 2x - c 2 u\ x + 2c 2 u x u 3x + J \3c 1 + c 2 )u x dx. (43) 

Removing the "obstruction" (the last term) requires c 2 = — 3c\. Substituting a solution, 
c i = |)C2 = —1, into (27) gives in (33). Substituting that solution into (42) and (43) 
yields 

/ = u 3 u x - 2u\ - 2uu x u 2x + u 2 u 3x - 2u x u ix (44) 
which is exact, together with its integral 

F = -u A — 2uu 2 x + u 2 u 2x + u\ x — 2u x u 3x , (45) 

which is in (33). Currently, CAS like Mathematical Maple 5 , and REDUCE cannot 
compute (43) as a sum of terms that can be integrated out and the obstruction. 

Example: Consider the following example in 2D 

/ = U X Vy - U 2x Vy - UyV X + U X yV X , (46) 

where u(x,y) and v(x,y). It is easy to verify by hand (via integration by parts) that / is 
exact. Indeed, / = DivF with 

F = (uv y — u x v y , —uv x + u x v x ). (47) 

3 We do not use integrable to avoid confusion with complete integrability from soliton theory [1, 2]. 
4 Variable t is parameter in subsequent computations. For brevity, we write u(x) instead of u(x;t). 
5 Future versions of Maple will be able to "partially" integrate such expressions [10, 19]. 
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As far as we know, the leading CAS currently lack tools to compute F. Two questions arise: 

(i) How can we tell whether / is the divergence of some differential function F? 

(ii) How do we compute F automatically and without multivariate integration by parts of 
expressions involving arbitrary functions? 

To answer these questions we use the following tools from the calculus of variations: the 
variational derivative (Euler operator), its generalizations (higher Euler operators or Lie- 
Euler operators), and the homotopy operator. 

Definition: The variational derivative (zeroth Euler operator), C^^, is defined [26, p. 246] 

by 



where the sum is over all the unordered multi-indices J [26, p. 95]. For example, in the 
2D case the multi-indices corresponding to second-order derivatives can be identified with 
{2x,2y,2z,xy,xz,yz}. Obviously, (-D) 2x = (-B x ) 2 = D*, {—D) xy = (_D X )(-D„) = D X D„, 
etc.. For notational details see [26, p. 95, p. 108, p. 246]. 

With applications in mind, we give calculus-based formulas for the variational derivatives 
in ID and 2D. The formula for 3D is analogous and can be found in [19]. 

Example: In ID, where x = x, for scalar component u(x) we have 6 

4(1) = Et-D*)*/- = #- - D x + D^ - Dl^- + ■ ■ ■ , (49) 
1 ' t^o du kx du du x du 2x du 3x 

In 2D where x = (x,y), we have for component u(x,y) 

£ (o,o) _^ V ( p ) k *( p ) k * 9 - d T) 9 D 9 

U{x > v) ~hh du kxXkyy -du du x du y 

+ +^ X B 9 + DJ-^ - D^^ , (50) 

ou 2x ou xy y ou 2y ou 3x 

Note that u kxXkyV stands for u xx ... xyy ... y where x is repeated k x times and y is repeated k y 
times. Similar formulas hold for components v,w, etc.. The first question is then answered 
by the following theorem [26, p. 248]. 

Theorem: A necessary and sufficient condition for a function / to be a divergence, i.e. there 
exists a differential function F so that / = DivF, is that 4°(L)(/) = 0- I 11 wor ds: the Euler 
operator annihilates divergences. 

If, for example, u(x) = (u(x),v(x)) then both C^) x Jf) and C^ljf) must vanish identically. 
For the ID case, the theorem says that a differential function / is exact, i.e. there exists a 
differential function F so that / = D X F, if and only if C^ x ^(f) = 0. 

6 Variable t in u(x; t) is suppressed because it is a parameter in the Euler operators. 
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Example: Avoiding integration by parts, we determine C\ and c 2 so that / in (42) will be 
exact. Since / is of order 4 in the (one) dependent variable u(x), the zeroth Euler operator 
(49) terminates at k — 4. Using nothing but differentiations, we readily compute 




(51) 



Note that the terms in u 2 u X) uu 3x , and u$ x dropped out. Hence, £> u ( x )(f) = leads to 
3ci + c 2 = 0. Substituting ci = |,c 2 = — 1 into (42) gives (44) as claimed. 

Example: As an example in 2D, we readily verify that / = u x v y — U2 X v y — u y v x + u xy v x from 
(46) is a divergence. Applying (50) to / for each component of u(x,y) = (u(x,y),v(x,y)) 



6.2 Removing Divergences and Divergence-equivalent Terms 

It is of paramount importance that densities are free of divergences for they could be moved 
into the flux J of the conservation law D t p + Div J = 0. We show how the Euler operator can 
be used to remove divergences and divergence-equivalent terms in densities. An algorithm 
to do so is given in [19]. The following examples illustrate the concept behind the algorithm. 

Definition: Two scalar differential functions, and f( 2 \ are divergence-equivalent if and 
only if they differ by the divergence of some vector V, i.e. ~ if and only if f^ — f^ = 
DivV. If a scalar expression is divergence-equivalent to zero then it is a divergence. 

Example: Any conserved density of the KdV equation (3) is uniform in rank with respect 
to W(u)—2 and W(d/dx) — l. The list of all terms of, say, rank 6 is 72. = {u 3 , u x , UU2 X , u Ax }. 
Now, terms = uu 2x and f < - 2 ' ) = —u 2 x are divergence-equivalent because — f < - 2 ' ) = 
u 2 x + uu 2x = D x (uu x ). Using (49), note that C^ x) (-u 2 x ) = 2u 2x and C^ x) (uu 2x ) = 2u 2x 
are equal (therefore, linearly dependent). So, we discard UU2 X . Moreover, u Ax = D x (m 3x ) is 
a divergence and, as expected, C^Ju^) = 0. So, we can discard u 4x . Hence, 7Z can be 
replaced by S = {u 3 ,u 2 } which is free of divergences and divergence-equivalent terms. 
Example: The list of non-constant terms of rank 6 for the Boussinesq equation (14) is 



Using (49), for every term U in U we compute V; = C ( u ^(ti) = (C u ( x )(ti), £ v ( x )(ti))- If 
Vj = (0, 0) then ti is discarded and so is v*. If Vj ^ (0, 0) we verify whether or not Vj is 



separately, we straightforwardly verify that £ u ( xy ){f) = and C> v ( X y){f) = 0- 



TZ = {P 2 u, f3u 2 , m 3 , v 2 , u x v, u 2 x , (3v x , uv x , (3u 2x , uu 2x , v 3x , u Ax }. 



(52) 
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linearly independent of the non-zero vectors v,-, j = 1, 2, • • • , % — 1. If independent, the term 
tt is kept, otherwise, ti is discarded and so is Vj. 

Upon application of (49), the first six terms in 1Z lead to linearly independent vectors vi 
through V6- Therefore, ti through te are kept (and so are the corresponding vectors). For 
t 7 = (3v x we compute v 7 = C^^flvx) = (0,0). So, t 7 is discarded and so is v 7 . For t 8 = uv x 

we get V8 = Cty x \(uv x ) = (v x , —u x ) = — v 5 . So, t$ is discarded and so is v 8 . 

Proceeding in a similar fashion, t 9 ,t 10 ,tn and t 12 are discarded. Thus, 1Z is replaced by 

S = {{3 2 u, (3u 2 , u 3 , v 2 , u x v , u x }, (53) 

which is free of divergences and divergence-equivalent terms. 

Example: In the 2D case, = (u x — u 2x )v y and = (u y — u xy )v x are divergence- 
equivalent since — = u x v y — u 2x v y — u y v x + u xy v x = Div (uv y — u x v y , —uv x + u x v x ). 
Using (50), note that C { °} x = ^ul, y) (f {2) ) = (-v xy - v xxy , -u xy + u xxy ). 

6.3 Higher Euler Operators 

To compute F = Div _1 (/) or, in the ID case, F = D~ 1 (f) = J f dx, we need higher-order 
versions of the variational derivative, called higher Euler operators [23, 26] or Lie-Euler 
operators [6]. The general formulas in terms of differential forms are given in [26, p. 367]. 
We give calculus-based formulas for the ID and 2D cases (see [19] for the 3D case). 

Definition: The higher Euler operators in ID (with variable x) 1 are 



oo 



4Vgy<-0^, (54) 

where (^j is the binomial coefficient. Note that the higher Euler operator for i — matches 
the variational derivative in (49). 

Example: The first three higher Euler operators in ID for component u(x) are 

£ (1) - — -2D — + 3D 2 — -4D 3 — + ••• 

u{x) ~ du x ZUx du 2x + 6Ux du 3x ^ x du 4x + ' 

4(1) = 3D,^ + 6D 2 -^ - 10D^ + ■ ■ ■ , (55) 

u(x) du 2x du 3x du 4x du 5x 

£ (3) = _iL_ 4D;c ^L + 10D 2 — -20D 3 ^- + ••• 
M{:r) du 3x x du 4x x du 5x x du 6x 

Definition: The higher Euler operators in 2D (with variables x, y) are 



oo oo 



4&j) = E E (•)(•) (—'D x ) kx ~ ix (—Dy) ky ~ iv - — . (56) 

k x =i x ky=i y \ l x ) \ l y ) V u k x xkyy 



7 Variable t in u(x; t) is suppressed because it is a parameter in the higher Euler operators. 
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Note that the higher Euler operator for i x =i y = matches the variational derivative in (50). 
Example: The first higher Euler operators for component u(x,y) are 

-(1,0) _ d ^ d d 2 d _ „ d 



u(x ' y > du x du 2x ou xy du 3x ou 2xy 

<S. = £ " 2D «^ " + 3D »^ + 2D * D ^ - ' • (57) 

- 19 2D - 2D - I 3D 2 - I ID D - I 

UKX ' y> du xy du 2xy ou x2y du 3xy ou 2x2y 



'«(*.!/) r7// 2 ,., ; X 9u 3sw v du 2x2y x du 4xy !l On 

The higher Euler operators are useful in their own right as the following theorem [23] in ID 
shows. 

Theorem: A necessary and sufficient condition for a function / to be an r th order derivative, 
i.e. there exists a scalar differential function F so that / = D X F, is that C^ x )(f) = for 
i = 0,l,...,r-l. 



6.4 Homotopy Operators 

Armed with the higher Euler operators we now turn to the homotopy operator which will 
allow us to reduce the computation of F = Div~ 1 (/), or in the ID case F = D~ 1 (/) = J f dx, 
to a single integral with respect to an auxiliary variable denoted by A (not to be confused 
with A in Section 4). Amazingly, the homotopy operator circumvents integration by parts 
(in multi-dimensions involving arbitrary functions) and reduces the inversion of the total 
divergence operator, Div, to a problem of single-variable calculus. 

As mentioned in Section 6.1, Div -1 is defined up to a divergence-free (curl) term. In 
3D, Div -1 is an equivalence class Div _1 (/) = F + V x K where K is an arbitrary vector 
differential function. The homotopy operator computes a particular K. 

The homotopy operator in terms of differential forms is given in [26, p. 372]. Below we 
give calculus-based formulas for the homotopy operators in ID and 2D which are easy to 
implement in CAS. The explicit formulas in 3D are analogous (see [19]). 

Definition: The homotopy operator in ID (with variable x) 8 is 

r l N d\ 

««(*)(/)= / £M/)[Au]-, (58) 
Jo j=1 A 

where Uj is the jth component of u and the integrand I Uj (f) is given by 

oo 

M/) = E D U«; <$(/))■ ( 59 ) 

i=0 

8 Variable t in u(x; t) is suppressed because it is a parameter in the homotopy operators. 
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The integrand involves the ID higher Euler operators in (54). In (58), iV is the number of 
dependent variables and I Uj (f)[Xu] means that in I Uj (f) we replace u(x) — > Au(:r), u x (x) — > 
\u x (x) : etc.. In practice, we first add the I Uj (f) and then scale the variables. The sum in 
Iuj{f) terminates at % — p — 1 where p is the order of Uj in /. 

Given an exact function /, the second question, that is how to compute F = D^ 1 (/) = 
/ fdx, is then answered by the following theorem [26, p. 372]. 

Theorem: For an exact function /, one has F = Ti. u (x)(f)- 

Thus, in the ID case, applying the homotopy operator (58) allows one to bypass integration 
by parts. As an experiment, one can start from some function F, compute / = D X F, then 
compute F = H u ( x )(f), and finally verify that F — F is a constant. 

Example: Using the homotopy operator (58) with (59), we recompute (45). Since / in (42) 
is of order p = 4 in Ui(x) = u(x), the sum in (59) terminates at i — p — 1 = 3. Hence, 

W) = <£)(/) + D x («4(l)(/)) + «^)(/)) + ^ (uC% } 



,K _ 2uD x pU + 3uB 2 pU - 4«D 3 pU + D x (u ° f - 3uD x (M- 
' du x x \duo x x \duzJ x \du Ax j x \ du 2x x \du a 



\du4xj J x \ du 3x \du4 X J J x \ du± x ) 

= -u 4 + 8u 2 U2x + Quu^x + Auul — 4D X (2u 2 u x + 3uu 3x ^j + D^ (u 3 + 8uu 2x ^j — 2D 3 (uu x ) 
= u 4 — 6uu 2 + 3u 2 u 2x + 2u\ x — Au x u 3x . (60) 

Formula (58), with u = Ui — u, requires an integration with respect to A : 

F = n u(x) (f) = / o 4(/)My 

= J (A 3 m 4 - 6X 2 uu 2 x + 3X 2 u 2 u 2x + 2\u\ x - A\u x u 3x ^j d\ 

= -u 4 — 2uu x + u 2 u 2x + u 2x — 2u x u 3x . (61) 

The crux of the homotopy operator method is that the integration by parts of a differential 
expression like (44), which involves an arbitrary function u(x) and its derivatives, can be 
reduced to a standard integration of a polynomial in A. 

Example: For a system with N = 2 components, u(x) = (ui(x),u 2 (x)) = (u(x),v(x)), the 
homotopy operator formulas are 

n u(x) (f) = £ (/„(/) + /„(/)) [Au] ^, (62) 

with 

oo oo 

W) = E°x(«4w ) (/)) and W) = E°x(«4w 1) (/))- ( 63 ) 

i=0 i=0 
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Example: Consider / = 3u x v 2 sin u — u x sin u — 6vv x cos u + 2u x u 2x cos u + Sv x v 2x , which is 
no longer polynomial in the two (N = 2) dependent variables u(x) and v(x). Applying Euler 
operator (48) for each component of u(x) = (u(x),v(x)) separately, we quickly verify that 
£u(x)(f) = and £$ x )(f) = 0. Hence, / is exact. Integration by parts (by hand) gives 

F — J f dx = Av 2 + u 2 cos-u — 3v 2 cos-u. (64) 

Currently, CAS like Mathematica, Maple 9 and Reduce fail this integration due to the presence 
of trigonometric functions. 

We now recompute (64) with the homotopy operator formulas (62) and (63). First, 

w) = <!,(/) +Dl «»,(/)) 

= 3uv 2 sin u — uu 2 x sin u + 2-u 2 cos u. (65) 

Second, 

9/ nn /9/\, n /9/^ 



c% V^x/ V ^2xy 

.2„„„„. , o„.2 



= — 6v cos-u + 8i; x . (66) 
Finally, using (62), 



dA 
T 



F=n u(x) (f)= [\i u (f) + i v (f))[\u] 

J 

= / (3X 2 uv 2 sin(A-u) — \ 2 uu 2 x sin(A-u)+2A-u 2 cos(A-u) — 6Ai> 2 cos(Aw)+8A^dA 

v 

= -u 2 cos-u — 3f 2 cos-u + Av 2 . (67) 

In contrast to the integration in (64) which involved arbitrary functions u(x) and v(x), in 
(67) we integrated by parts with respect to the auxiliary variable A. Any CAS can do that! 

We now turn to inverting the Div operator using the homotopy operator. 

Definition: We define the homotopy operator in 2D (with variables x, y) through its two 

components (H^x y)(f) Wu(x y)(f)) ■ The re-component of the operator is given by 



with 



E^ } (/)[Au]^, (68) 

JO j=1 A 
OO OO / -I I • \ 

m) = EE rrfrr D * D * ^<<^) ■ (69) 



i x =0 j u =0 



9 Version 9.5 of Maple can integrate such expressions as a result of our interactions with the developers. 
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Likewise, the y-component is given by 



H%Jf) =J E^(/)[Au]f, (70) 



with 

4f (/) = E E ( i^fr-l D " D ^ (^-^S^/)) • ( 71 ) 

Integrands (69) and (71) involve the 2D higher Euler operators in (56). The question how 
to compute F = F 2 ) = Div~ (/) is then answered by the following theorem. 

Theorem: If / is a divergence, then F = (F U F 2 ) = Div"^/) = (n%^f),H^ y) (f)). 
The superscripts (x) and (y) remind us which components of F we are computing. As a 
matter of testing, we can start from some vector F and compute / = DivF. Next, we 
compute F = (F 1 ,F 2 ) = (n% ty) (f),n% ty) (f)) and, finally, verify that K = F - F is 
divergence free. 

Example: Using (46), we show how application of the 2D homotopy operator leads to (47), 
up to a divergence free vector. Consider / = u x v y — u 2x v y — u y v x + u xy v x , which is the 
divergence of F in (47). In order to compute Div -1 (/), we use (69) for the u component in 
u(x,y) = (u(x,y),v(x,y)) : 

i { :\f) = ^(/i+D.Kiaij + ^^^y/)) 

U (du x 2 ^ >x du 2x ^ y du X y) + ^ x { U du 2x ) + 2 Dy ( U du xy ) 
11 

= UVy + ~UyV x - U x Vy + ~UV X y . (72) 

Similarly, for the v component of u(x,y) = (u(x,y),v(x,y)) we get 

4 x) (f) = = « J£ = -u y v + u xy v. (73) 

Hence, using (68), 

*i = = [ {4 x) (f) + 4 x) (f)) [Au] 

= X (iLVy + ^UyV X - U X Vy + ^UV X y ~ UyV + U X yV^ (IX 

= ^UVy+-UyV X --U X Vy + -UV X y--U y V+-U X yV. (74) 

Without showing the details, using (70) and (71) we compute analogously 
F2 = H^Jf) = [ (#>(/) + #>(/)) [Au] ^ 

19 



J x 



dX 

T 



= A (-uv x - ^uv 2x + -jU x v x + u x v - u 2x vj d\ 

11 1 1 1 

= ~-uv x - -uv 2x + -u x v x + -u x v - -u 2x v. (75) 

Now we can readily verify that the resulting vector 

F = l Fl \ = / \uVy + \ U yV X - \U X Vy + \uV X y ~ \UyV + ^yV \ ^ 

\ F 2 ) \ -\uv x - \uv 2x + \u x v x + \u x v - \u 2x v J 
differs from F = (uv y — u x v y , —uv x + u x v x ) by a divergence-free vector 

K = F _ f= | ^ 1 i = ( 2 uv y ~ ~i u y Vx ~ 2 UxV v ~ i UVx v = 2 u v v ~ 2 u xy v \ 
\K 2 J \ -\uv x + \uv 2x + \u x v x - \u x v + \u 2x v )' 

As mentioned in Section 6.1, K can be written as (D^, —D x (p) with <fi = ^uv — \uv x — \u x v. 



7 Computation of Conservation Laws 

In this section we apply the Euler and homotopy operators to compute density-flux pairs 
for the four PDEs in Section 3. Using the KdV equation (3) we illustrate the steps for a 
ID case (space variable x) and dependent variable u(x,t). Part of the computations for this 
example were done in the previous sections. The Boussinesq equation (14), also in ID, is 
used to illustrate the case with two dependent variables u(x,t) and v(x, t) and an auxiliary 
parameter with weight. 

The LL equation (7), still in ID, has three dependent variables u(x, t),v(x, t) and w(x, t). 
The three coupling constants have weights which makes the computations more cumbersome. 
However, the multi-uniformity of (7) helps to reduce the complexity. The shallow-water wave 
equations (10) illustrate the computations in 2D (two space variables) and four dependent 
variables. Again, we can take advantage of the multi-uniformity of (10) to control "expression 
swell" of the computations. 

We use a direct approach to compute conservation laws, D t p + Div J = 0, of polynomial 
systems of nonlinear PDEs. First, we build the candidate density p as a linear combination 
(with constant coefficients q) of terms that are uniform in rank with respect to the scaling 
symmetry of the PDE. In doing so, we dynamically remove divergences and divergence- 
equivalent terms to get the shortest possible candidate density. 

Second, we evaluate D t p on solutions of the PDE, thus removing all time derivatives 
from the problem. The resulting expression, E = — D t p, must be a divergence of the as yet 
unknown flux. Thus, we compute £^L(F) and set the coefficients of like terms to zero. 
This leads to a linear system for the undetermined coefficients q. If the given PDE has 
arbitrary constant parameters, then the linear system is parametrized. Careful analysis of 
the eliminant is needed to find all solution branches, and, when applicable, the conditions 
on the parameters. For each branch, the solution of the linear system is substituted into p 
and E. 
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Third, we use the homotopy operator H u ( x ) to compute J = Div~ 1 (£'). The computations 
are done with our Mathematica packages [18]. Recall that J is only defined up a curl term. 
Removing the curl term in J may lead to a shorter flux. Inversion of Div via the homotopy 
operator therefore does not guarantee the shortest flux. 



7.1 Conservation Laws for the KdV Equation 

In (33) we gave the first three density-flux pairs of (3). As an example, we will compute p^ 
and J (3) . The weights are W(d/dx) = 1 and W(u) = 2. Hence, density p*- 3 -* has rank 6 and 
flux J® has rank 8. The algorithm has three steps: 
Step 1: Construct the form of the density 

Start from V = {«}, i.e. the list of dependent variables (and parameters with weight, if 
applicable). Construct the set M. of all non-constant monomials of (selected) rank 6 or less 
(without derivatives). Thus, M. = {u 3 , u 2 , u}. Next, for each monomial in A4, introduce the 
needed ^-derivatives so that each term has rank 6. Since W{d/dx) = 1, use 

2 u 2 O^u 

-^- = {2uu x ) x = 2u 2 x + 2uu 2xi g^A =Uix - ( 78 ) 

Ignore the highest-order terms (typically the last terms) in each of the right hand sides of 
(78). Augment M. with the remaining terms (deleting numerical factors) to get 1Z = {u 3 , u 2 x }. 

Here S — 1Z since 1Z is free of divergences and divergent-equivalent terms. Linearly 
combine the terms in S with constant coefficients to get the candidate density: 

p = civ 3 + c 2 ul, (79) 

which is (27). 

Step 2: Compute the undetermined constants q 

Compute D t p and use (3) to eliminate u t and u tx . As shown in (28), this gives 

E = — D t p = 3ciu 3 u x + 2>c\u 2 u^ x + 2c2« 3 + 2c2Uu x U2 X + 2c2U x u^ x . (80) 

Since E = —D t p = D X J, the expression E must be exact. As shown in (51), C^) X JE) = 
— 6(3ci + C2)u x u 2x = leads to ci = |,C2 = — 1 and upon substitution into (79) to 

P = \u*~ul (81) 

which is in (33). 

Step 3: Compute the flux J 

Substitute Ci = |,c 2 = — 1 into (80) to get 

E = u 3 u x - 2w 3 - 2uu x u 2x + u 2 u 3x - 2u x u Ax . (82) 
As shown in (61), application of (58) with (59) to (82) gives 

J = -M 4 — 2uu 2 x + u 2 u 2x + u\ x — 2u x u 3x , (83) 

which is J( 3 ) in (33). 
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7.2 Conservation Laws for the Boussinesq Equation 

The first few density-flux pairs of (14) were given in (35). Equation (14) has weights W(-^) = 
1, W(u) — W(/3) = 2, and W(v) = 3. We show the computation of p^ and of ranks 6 and 
7, respectively. The presence of the auxiliary parameter j3 with weight complicates matters. 
Step 1: Construct the form of the density 

Augment the list of dependent variables with the parameter (3 (with non-zero weight). Hence, 
V = {u,v,/3}. Construct M. = {/3 2 u, /3u 2 , flu, (3v,u 3 ,u 2 ,u,v 2 ,v,uv}, which contains all non- 
constant monomials of (chosen) rank 6 or less (without derivatives). Next, for each term in 
M., introduce the right number of ^-derivatives so that each term has rank 6. For example, 

d 2 f3u d 2 u 2 2 d 4 u d(uv) 

-^- = pu 2x , -^ = 2u x + 2uu 2x , t^=u 4x , — = u x v + uv x , etc.. (84) 

Ignore the highest-order terms (typically the last terms) in each of the right hand sides 
of (84). Augment M. with the remaining terms, after deleting numerical factors, to get 
1Z = {(3 2 u, flu 2 , u 3 , v 2 , u x v, u x , (5v x , uv x , (3u 2x , uu 2x , v 3x , u 4x } as in (52). As shown in (53), 
removal of divergences and divergence-free terms in TZ leads to S = {{3 2 u, j3u 2 , u 3 , v 2 , u x v, u 2 .}. 
Linearly combine the terms in S to get 

p = c\fl 2 u + c 2 /3u 2 + c 3 -u 3 + c 4 v 2 + c 5 u x v + c 6 u 2 x . (85) 

Step 2: Compute the undetermined constants q 

Compute E = — D t p and use (14) to eliminate u t ,u tx , and v t . As shown in (32), this gives 
E = (ci/3 2 + 2c 2 (3u + 3c 3 u 2 )v x + (c 5 v + 2c 6 u x )v 2x + (2c 4 v + c b u x )((3u x - 3uu x - au 3x ). (86) 

E = —D t p = D x J must be exact. Thus, apply (48) and require that C^) X JE) = C^) X \(E) = 0. 
Group like terms. Set their coefficients equal to zero to obtain the parametrized system 

j3(c 2 — C4) =0, C3 + C4 = 0, C5 = 0, ac$ = 0, /^cs = 0, ac 4 — c% = 0. (87) 

Investigate the eliminant of the system. Set c± — 1, to obtain the solution 

d = 1, c 2 = c 4 , c 3 = -c 4 , c 5 = 0, c 6 = ac 4 , (88) 

which is valid irrespective of the values of a and (3. Substitute (88) into (85) to get 

p = [3 2 u + c 4 {(3u 2 -u 3 + v 2 + au 2 x ). (89) 

The density must be split into independent pieces. Indeed, since c 4 is arbitrary, set c 4 = 
or c 4 = 1, thus splitting (89) into two independent densities, p = fi 2 u and 

p = f3u 2 - u 3 + v 2 + au 2 x , (90) 

which are p^ and p^ in (35). 
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Step 3: Compute the flux J 

Compute the flux corresponding to p in (90). Substitute (88) into (86). Take the terms in 
C4 and set C4 = 1. Thus, 

E = 2(3uv x + 2(3u x v — 3u 2 v x — 6uu x v + 2au x V2 X — 2au 3x v. (91) 

Apply (62) and (63) to (91). For E of order 3 in u(x), compute 

Iu{E) = uC^E) + D x (uC%{E)) + Dl (uC® x) (E)) 

dE n ^ ( dE \ o - / dE \ ^ f dE n ^ / dE \\ - / dE \ 
= u- 2uD x + 3uD 2 +D X [u- 3uD x + Y)l [u- 



du x " \du 2x x \du 3x l \ 9u 2x \du 3x JJ x \ du 3xJ 



= 2(3uv — 6u 2 v — Aauv 2x + 6aD x (uv x ) — 2aD x (uv) 

= 2j3uv — 6u 2 v — 2au 2x v + 2au x v x . (92) 
With E is of order 2 in v{x), subsequently compute 

LAB) = ^(E) + D. (*C» ){ B» = ,f - 2.D, (g + D, („g 

= 2/3-ut> — 3-u 2 f — Aauv 2x + 2aD x (u x v) = 2(3uv — 3u 2 v — 2au 2x v + 2au x v x . (93) 
Formula (62) with u(x) = (u(x),v(x)) requires an integration with respect to A. Hence, 

J = H U{X) {E) = j\lu(E) + I v (E))[\u] ^ 

= J (Af3\uv — 9\ 2 u 2 v — Aa\u 2x v + Aa\u x v x ^j dX 

= 2j3uv — 3u 2 v — 2au 2x v + 2au x v x , (94) 

which is j( 4 ) in (35). Set f3 — 1 at the end of the computations. 

7.3 Conservation Laws for the Landau-Lifshitz Equation 

Without showing full details, we will compute the six density-flux pairs (36) for (7). The 
weights for (7) are W(u) = W(v) = W(w) = a, W(a) = W(0) = W(i) = 2, W(§- t ) = a + 2, 
where a is an arbitrary non-negative integer or rational number. This example involves three 
constants (a, (3, and 7) with weights which makes the computations quite cumbersome unless 
we take advantage of the multi-uniformity and the cyclic nature of (7). 
Step 1: Construct the form of the density 

Select a small value for a, for example a — |, and compute a density of, say, rank |. Start 
from V = {u, v, w, a, /3, 7}, i.e. the list of dependent variables and the parameters with 
weight. Construct the set A4 of all non-constant monomials of the (selected) rank | or less 
(without derivatives). Thus, M. = {u, v, w}. Next, for each of the monomials in A4, introduce 
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the appropriate number of ^-derivatives so that each term has rank |. No ^-derivatives are 
needed and the removal of divergences or divergence-equivalent terms is irrelevant. Linearly 
combine the terms in TZ—S — {u, v, w} to get a candidate density: 

p — ciu + c 2 v + c 3 w. (95) 

It suffices to continue with p = C\U. The remaining terms follow by cyclic permutation. 
Step 2: Compute the undetermined constants q 

Compute E = —D t p = —D t (ciu) = —c\u t . Use (7) to eliminate u t . Hence, 

E = -c 1 (vw 2x - wv 2x + (7 - P)vwj, (96) 



which is the opposite of the right hand side of the first equation in (7). Since E— — D t p = D x J, 

■j,(x) ( 



the expression E must be exact. Obviously C^K(E) = 0. Compute 



Set the latter expressions identically equal to zero. Solve C\(f3 — 7) = for c\ 7^ 0. Set C\ — 1 
to get p — u, subject to the condition f3 — 7, which confirms the result in (36). 

Step 3: Compute the flux J 

Substitute c± — 1 and (5 — 7 into (96) to get i? = w%-fW2i. Apply the homotopy operator 
(58) with (59) to E. In this example u(x) = (u(x),v(x),w(x)). Obviously, I U (E) = 0, since 
there is no explicit occurrence of u. Compute 

I v (E)=vC i ^ ) (E)+D x (vC%(E))=v— -2vD x \—j+D x \v—j=v x w-vw x , (98) 

and, analogously, I w (E) = v x w — vw x . Finally, compute 

dX 



J = H U(X) (E) = j\lu(E) + I V (E) + I w (E))[Xu] 



I 



A 

2A (v x w — vw x ) dX = v x w — vw x , (99) 



which is in (36). 

To compute density p*- 4 -* which is quadratic in u, v, and w, we would start from rank | and 
repeat the three steps. Step 1 would lead to p = c\u 2 + c 2 uv + c 3 v 2 + c^uw + c 5 vw + c G w 2 . 
Steps 2 and 3 would result in c\ = c 3 = c 6 = 1. So, p = u 2 + v 2 + w 2 and J = 0, which agrees 
with and j( 4 ) = in (36). 

Starting with rank 1, Step 1 would generate a polynomial density with the 15 terms that are 
homogeneous of degree 4. Steps 2 and 3 would result in p^ and = in (36). 
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The computation of and is cumbersome and long. We only indicate the strategy and 
give partial results. Notice that p^ would have rank | if a — | and the candidate density 
would have, amongst others, all homogeneous terms of degree nine, which is undesirable. It 
turns out that a = 1 is a better choice to compute p*- 6 -* , which then has rank 4. 

Start from V as above. Construct the set M. = {u 4 , u 3 , u 2 , u, t> 4 , v 3 , • • • , u 3 v , u 3 w, • • • , 
u 2 v 2 , u 2 w 2 , • • • , u 2 vw, v 2 uw, • • • , au 2 , f3u 2 , • • • , auv, • • • , •yvw, • • • , uvw, uv, • • • , •yw} of all non- 
constant monomials of rank 4 or less (without derivatives). Next, for each of the 61 mono- 
mials in Ai, introduce the needed ^-derivatives to make each term rank 4. Use, for example, 

d 2 u 2 ^ 2 du 3 n 2 d 4 u du 2 v 2 9au , inn . 

-Q^T = 2u x + 2uu 2x, -r^ = 3uu x , — = u 4x , -^— = 2uu x v + uv x , -^— = au x . (100) 

Ignore the highest-order terms (typically the last terms). Augment M. with the remaining 
terms (deleting numerical factors). Next, remove all divergences and divergent-equivalent 
terms to get S with 47 terms (not shown). Linearly combine the terms in S to get the 
candidate density: 

p = CiCm + C 2 [3ll + C37M + C4M + C5Q;Mf + Cq(3uV + Cj^/UV + CsM f + • • • + C12M f 

H h c 2 2U 2 vw + c 23 uv 2 w H h c 25 cra; 3 + c 2e /3w 3 + c 27 yw 3 H h c 29 uvw 2 

+C 30 V 2 W 2 H h C33W 4 + C 34 MMx^ + c 35 u x v 2 + c 36 uu x w + c 37 u x vw + c 38 u x w 2 + c 39 u 2 x 

+c 40 uv x w + c 41 vv x w + c 42 v x w 2 + c 43 u x v x + c u vl + c 45 u x w x + c m v x w x + c i7 wl, (101) 

where only the most indicative terms are explicitly shown. 

Compute D t p, use (7) to eliminate u t , v t , w t , u tx , v tx , and w tx , and require that C^ X ^(E) = 

Cf^(E) =£wlx)(E) = which leads to a linear system of 121 equations (not shown) for c\ 
through C47. Substitute the solution (not shown) into (101) to obtain 

P = («c 2 5 + (3c 26 + 7c 27 )(m 2 + v 2 + w 2 ) + ^c 30 (w 2 + v 2 + w 2 ) 2 

+c 47 (ul + v 2 x + w 2 x + (7 - a)u 2 + (7 - (3)v 2 ) , (102) 
where c 25 , c 26 , c 27 , c 30 , and C47 are arbitrary. Split the density into independent pieces: 

2 , 2 1 2 / 2 1 2 1 2\2 

P = < + v 2 x + w 2 x + (7 - a)u 2 + (7 - (3)v 2 , (103) 

which are p^\ p^\ and p^ 6 ^ in (36). Use the homotopy operator to compute the flux corre- 
sponding to (102): 

J = 2c A1 [{vw x - wv x )u 2x + (wu x - uw x )v 2x + (uv x - vu x )w 2x 

+ ((3 — y)vwu x + (7 — a)uwv x + (a — (3)uvw x ^j. (104) 

J(4) = j(5) = o 

since the terms in C25, c 2 6, c 2 7, and C30 all dropped out. Set C47 = 1 to get 

J( 6 ) in (36). 
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7.4 Conservation Laws for the Shallow Water Wave Equations 

The SWW equations (10) admit weights (19) in which W(h) and W(Q) are free. The fact 
that (10) is multi-uniform is advantageous. Indeed, we can construct a candidate p which is 
uniform in rank for one set of weights and, subsequently, use other choices of weights to split 
p into smaller pieces. This "divide and conquer" strategy drastically reduces the complexity 
of the computations as was shown in [19]. 

The first few densities and fluxes were given in (41). We compute p^ and J^. Note 
that p^ has rank 3 if we select W(h) = a = 1 and W{Vt) = b = 2. However, p^ has rank 4 
if we take W(h) = a = and W(Q) = b = 2. If we set W(h) = 0, and W(Q) = 3 then p^ has 
rank 7. So, the trick is to construct a candidate density which is scaling homogeneous for a 
particular (fixed) choice of a and b in (20) and split the density based on other choices of a 
and b. 

Step 1: Construct the form of the density 

Start from V = {u,v,9, h,Q}, i.e. the list of variables and parameters with weights. Use 
(20) with a — 1, b — 2, to get M — {flu, Qv, . . . , u 3 , v 3 , . . . , u 2 v, uv 2 , . . . , u 2 , v 2 , . . . , u, v , 9, h} 
which has 38 monomials of (chosen) rank 3 or less (without derivatives). 

All terms of rank 3 in M. remain untouched. To adjust the rank, differentiate each 
monomial of rank 2 in M. with respect to x ignoring the highest-order term. For example, 
in = 2uu x , the term can be ignored since it is a total derivative. The terms u x v and 
— uv x are divergence-equivalent since = u x v + uv x . Keep u x v. Likewise, differentiate 
each monomial of rank 2 in Ai with respect to y and ignore the highest-order term. 

Produce the remaining terms for rank 3 by differentiating the monomials of rank 1 in 
Ai with respect to x twice, or y twice, or once with respect to x and y. Again ignore 
the highest-order terms. Augment the set Ai with the derivative terms of rank 3 to get 
1Z = {flu, Qv, • • • , uv 2 , u x v, u x 8, u x h, ■ ■ ■ , u y v, u y 9, ■ ■ ■ , y h} which has 36 terms. 

Use the "divide and conquer" strategy [19] to select from 1Z the terms which are of 
ranks 4 and 7 using the choices a = 0, b = 2 and a = 0, b = 3, respectively. This gives S = 
{Q9,u x 9,u y 9,v x 9,v y 9}, which happens to be free of divergences and divergence-equivalent 
terms. So, no further reduction is needed. Linearly combine the terms in S to get 



p = ciVl9 + c 2 u x 9 + c 3 u y 9 + c 4 v x 9 + c 5 v y 9. 



(105) 



Step 2: Compute the undetermined constants q 

Compute E = —D t p and use (10) to remove all time derivatives: 




c<i9{uu x + vu y — 2Qv + \h9 x + 9h x ) x + c^9{uu x + vu y — 2Qv + \h9 x + 9h x ) y 
+ c±9{uv x + Wy + 2Qu + \h9 y + 9h y ) x + c 5 9(uv x + vv y + 2Qu + \h9 y + 9h y ) 
+ (ciQ + c 2 u x + QUy + c 4 v x + c 5 v y )(u9 x + v9 y ). 



(106) 
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Require that C^ y) {E) = C^ y) {E) =C™ y) (E) = C™ y) (E) = 0, where, for example, C^ y) 
is given in (50). Gather like terms. Equate their coefficients to zero to obtain 

ci + 2c 3 = 0, c 2 = c 5 = 0, ci-2c 4 = 0, c 3 + c 4 = 0. (107) 

Set Ci = 2. Substitute the solution 

ci = 2, c 2 = 0, c 3 = -1, c 4 = 1, c 5 = 0. (108) 

into p to obtain 

p = 2Q9 -u y 6 + v x 6, (109) 



which is p( 5 ) in (41). 



Step 3: Compute the flux J 

Compute the flux corresponding to (109). To do so, substitute (108) into (106) to obtain 

E = 6{U X V X - UyVy + UV 2X - U 2 yV - U X Uy + V x Vy ~ UU X y + VV X y 

+2Qu x + 2Qv y + \6 x h y — \6 y h x ) + 2Qu9 x + 2Qv9 y — uu y 9 x 

—UyvOy + uv x 6 x + vv x 9 y . (HO) 

Apply the 2D homotopy formulas in (68)-(71) to E = Div J = D x Ji + D y J 2 . So, compute 
IP(E) = uC^ y) (E) + D x (uC^ y) (E)) + ±D S («4(i)(^ 

(—-2D ( dE \ D f— ^ D f — ^ Id ( — ^ 

\0u x x V 9^2x7 V \ du xyJJ X \du 2 J 2 v \du xy ) 
= uv x 6 + 2Vtu6 + ^u 2 e y - uu y 6. (Ill) 



Similarly, compute 



/(*)(£) = w „0 + ^0,, + (H2) 

4 x) (s) = ^e 2 h y + 2n u e -uu y e + uv x e, (113) 

(£) = -±60 W A. (114) 



Next, compute 



■M«) = <L,y)(^) 



= f^^) + 4 X) (E) + I<T\E) + tf\Ej) [Au] f 

= J (aXVLuO + A 2 (suv x 9 + 7jU 2 6 y - 2uu y d + w„0 + ^v 2 9 y 

+ l -9 2 h y - \ee y hj)d\ 

2 11111 

= 2Qu9 uuy6 + uv x 6 + -WyO + -u 2 6 y + -v 2 9 y h99 y + -h y 9 2 . (115) 

3 3 6 6 6 6 
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Likewise, compute I^(E)J^(E),I^ y) (E) : and lj?\E). Finally, compute 
^(u) = H^l y) {E) 

= + li v \E) + I<?\E) + [Au] ^ 

= 2Qv9 + ^w x - vu v 9 - -uu x 9 - -u 2 9 x - -v 2 9 x + -h99 x - -h x 9 2 . (116) 
3 3 6 6 6 6 



Hence, 



Ji \ _ 1 / 12M - 4m* y + 6uv x 9 + 2w„0 + u 2 9y + v 2 ^ - /i00„ + h y 9 2 \ , . 
J 2 J ~ 6 1 12fiw^ + 4w x - 6to^ - 2«M - u 2 9 x - v 2 9 x + - /i^ 2 i ' ^ ' 



which matches J^ 5 -* in (41). 



8 Conclusions 

The variational derivative (zeroth Euler operator) provides a straightforward way to test 
exactness which is key in the computation of densities. The continuous homotopy operator 
is a powerful, algorithmic tool to compute fluxes explicitly. Indeed, the homotopy operator 
allows one to invert the total divergence operator by computing higher variational derivatives 
followed by a one-dimensional integration with respect to a single auxiliary parameter. The 
homotopy operator is a universal tool that can be applied to problems in which integration 
by parts (of arbitrary functions) in multi- variables is crucial. 

To reach a wider audience, we intentionally did not use differential forms and the abstract 
framework of variational bi-complexes and homological algebra. Instead, we extracted the 
Euler and homotopy operators from their abstract setting and presented them in the language 
of standard calculus, thereby making them widely applicable to computational problems in 
the sciences. Our calculus-based formulas can be readily implemented in CAS. 

Based on the concept of scaling invariance and using tools of the calculus of variations, 
we presented a three-step algorithm to symbolically compute polynomial conserved densities 
and fluxes of nonlinear polynomial systems of PDEs in multi-spatial dimensions. The steps 
are straightforward: build candidate densities as linear combinations (with undetermined 
constant coefficients) of terms that are homogeneous with respect to the scaling symmetry 
of the PDE. Subsequently, use the variational derivative to compute the coefficients, and, 
finally, use the homotopy operator to compute the flux. 

With our method one can search for conservation laws in chemistry, physics, and engi- 
neering. Symbolic packages are well suited to assist in the search which covers many fields 
of application including fluid mechanics, plasma physics, electro-dynamics, gas dynamics, 
elasticity, nonlinear optics and acoustics, electrical networks, chemical reactions, etc.. 
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